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Abstract Research within the field of multiscale modelling seeks, amongst 
other questions, to reconcile atomistic scale interactions with thermodynam¬ 
ical quantities (such as stress) on the continuum scale. The estimation of 
stress at a continuum point on the atomistic scale requires a pre-defined ker¬ 
nel function. This kernel function derives the stress at a continuum point 
by averaging the contribution from atoms within a region surrounding the 
continuum point. Commonly the kernel weight assignment is isotropic: an 
identical weight is assigned to atoms at the same spatial distance, which is 
tantamount to a local constant regression model. In this paper we employ a 
local linear regression model and leverage the mechanism of automatic kernel 
carpentry to allow for spatial averaging adaptive to the local distribution of 
atoms. As a result, different weights may be assigned to atoms at the same 
spatial distance. This is of interest for determining atomistic stress at stack¬ 
ing faults, interfaces or surfaces. It is shown in this study that for crystalline 
solids, although the local linear regression model performs elegantly, the ad¬ 
ditional computational costs are not justified compared to the local constant 
regression model. 
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1 Introduction 

Atomistic to continuum coupling seeks to link atomistic extensive quantities 
(for example mass, momentum, energy) to a continuum quantity of interest 
(for example, stress, heat flux). This link has previously been achieved by 
defining a space averaging volume (as given by a kernel function muisi 
I18] l over which the atomistic extensive quantities are smoothed to compute 
the desired continuum quantity. The kernel function adheres to the standard 
constraints of a proper smoothing function Section 4.2]. 

Recently, interest in modelling atomistic stress [MlH] has been growing 
primarily due to the rise in computational power and the development of 
atomistic to continuum multiscale methods nnziinn]- Furthermore, attempts 
have been made to determine the kernel function using statistical arguments 
[SIEQ], correlation functions from statistical mechanics [26] or data-driven 
methods from machine learning EaUH]. A different expression for the atom¬ 
istic stress, for the particular study of free surfaces, has been investigated [H 
124] based upon a force-area concept. In previous work, the influence of the 
kernel decays equally in all spatial directions with increasing distance from 
the centre of the kernel. In contrast, little is currently known about the effect 
of allowing the kernel weighting to vary in different dimensions according to 
the local distribution of the data. The possibility of a-priori defining such 
a kernel function according to a given anisotropic distribution of matter is 
discussed in Enio]. In this paper we address this gap in the literature by 
investigating a non-parametric kernel regression method m, known as local 
linear regression (LLR), that permits the weighting given by the kernel to 
adapt to the local distribution (or density) of the data. In contrast to [2T1 
120] our method is entirely data-driven, with the model constructed directly 
from the atomistic data itself. 

In this work we take the observation made by |28| as our starting point, 
namely that the standard definition of virial stress is equivalent to a local 
constant regression (LCR) model, and investigate the benefit of extending 
the LCR model to permit the kernel weighting to vary along different spatial 
directions. LCR approximates the regression function by a (local) constant. 
A well-known issue with LCR is that boundary problems can detract from 
the global performance of the estimator: at the boundary of the predictor 
space the kernel neighbourhood is asymmetric, which can cause a substan¬ 
tial increase in the bias of the estimate. For example, a kernel at a point at 
the boundary will be truncated in its extent - those neighbouring data-points 
within this truncated extent will disproportionately influence the resulting 
predicted average by pushing or pulling the average above or below the value 
it would ordinarily assume if there were no boundary at that point. This 
is also true of irregular or non-uniform interior regions. Both situations can 
frequently be found in crystalline solids: namely at stacking faults and the 
surfaces of flaws. This bias somewhat restricts application of the LCR esti¬ 
mator to values in the interior of the material. 

In this contribution we study the applicability of LLR to the modelling of 
atomistic stress, which is a known technique for reducing the bias significantly 
at the boundaries at a modest cost in variance. Rather than making a local 
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constant approximation, as for LCR, the LLR model instead fits a linear 
regression line through the observations in the neighbourhood of each target 
data-point. It can be shown that fitting a linear regression estimate at each 
target point is the same as computing the response by taking a weighted 
summation of surrounding data-points, weighted by an equivalent or effective 
kernel which has been automatically adapted to the density of the samples. 
This phenomenon is known as automatic kernel carpentry in the statistical 
literature [3 Section 6.1.1]. In this paper we are particularly interested in 
applying LLR to materials that have surfaces and flaws in their lattice where 
the data-adaptive property of the equivalent kernel should be particularly 
noticeable. Our main point of comparison will be the LCR model advocated 
by [55]. 

The organization of the paper is as follows. Section [2] in cludes a brief 
review of the virial and Hardy definitions of stress. Section]^ introduces the 
LLR regression model and describes the formulation of the equivalent kernel 
which we incorporate into the atomistic stress definition. The performance 
of these kernels is tested in Section]^ Finally Section [^discusses our findings 
and provides suggestions for future work in this area. 


2 Atomistic stress 

The Irving and Kirkwood wm formulation of stress relates continuum quan¬ 
tities such as stress and heat flux to microscopic kinematics and kinetics (for 
an elaborate discussion see wm)- Given their definition, we can derive the 
virial stress (t„ at a continuum point x by introducing a space averaging 
volume 17 (x) centred at point x and a uniform kernel tfi = 1/I7(x) that has 
the dimension of inverse volume: 

N ( 1 ^ 

(T„(x) = - Y.i=l V'iCTi = - ^ V'i "liUi G Ui -I- — ^ fy (g) Xy -I- . . . 

i=l \ ' i = 

where N is the number of atoms contained within the averaging volume, 
denotes the mass of atom z, and is the relative velocity of the given atom 
to the mean velocity of the N atoms. The position vector of atom i is given 
by Xi, with x.ij = x^ — Xj denoting the distance between atoms i and j, and 
tij specifying the force on atom i due to its pair interaction with atom j. 
The stress in Equation may be equated to the virial theorem [51I15LI16| by 
replacing 17 (x) with the total volume of the system. 

Based upon the Irving and Kirkwood procedure we can obtain the Hardy 
definition wm of stress by introducing the kernel function tfi = ip{\x — Xi\): 

(Thix) = - (^ipi-miU^ (g) Ui ^ Bijiij (g ^ij + .. .y (2) 

The kernel function ipi defines the extent of the space averaging volume 
surrounding a continuum point [201 Section 4.2]. The Hardy dehnition of 
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stress introduces a bond function ' 0 (|x — + Ax^ |) dA between 

atoms i and j. 


3 Local Linear Kernel Regression 


Kernel regression m models the regression relationship between an ex¬ 
planatory variable X and a response variable Y as: 

y, = m{X,) + (3) 

where i = 1... N and m(») is the regression function, N is the number of 
data points and et is independently and identically distributed zero mean 
noise. The goal of kernel regression is to estimate the mean response curve m 
in the regression relationship. If we assume m is smooth, namely observations 
close to X contain information about the value of m at X, then we can use 
the P-term Taylor expansion of m at X (if X is near to the sample Xi) to 
estimate the value of the function at X: 

m{X,) « m(X) +m'(X)(X, - X) + ^m"(X)(X, - X f + .... (4) 


The parameters /3fc, with k = 0 ... P, can be estimated by solving the follow¬ 
ing least squares optimisation problem: 

N 

rnin^ (X, - (3o - /3i(X, - X) - /32(X, - X)^ - ...)" X(|X, - X|), (5) 

Pfc . _ 

1 — 1 


where the kernel function K gives nearby samples higher weight than more 
distant samples. The kernel function is constrained to be non-negative, sym¬ 
metric and uni-modal. 

The previous treatment can be generalised into a multivariate (in this case 
3-dimensional) problem by considering a dataset of quadruples (xi,yi),..., 
{xn,Yn) with Xi = [xi,yi, Zi]"^. In this case the least squares optimisation 
problem (Equation]^ may be stated as: 


N 

min (Vi — m(x) — Vto(x) -(x^ — x) 

- (x* - x) • ^'Hm(x) -(x, - x) - ... ^ X(|x, - x|), ( 6 ) 


/32 


in which V(*) and 'H(») give the gradient and the Hessian, respectively. The 
symbol (•) gives the dot product and ((8)) the dyadic product 

We introduce Voigt’s notation, for example, given a symmetric 3x3- 
matrix A we represent it as a vector v^A} = [An, A 22 , A 33 , A 23 , A 13 , ^ 12 ]^. 
By doing so we can rewrite Equation rain matrix form as [zuaiin] 

's = argmin(Y — Bs)^W(Y — Bs) 

S 


( 7 ) 
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with Y = s = [^o,/3i^, vnj/Ja}^, • • W = diag[i4:(|xi - 

x|),i 4 r(|x 2 — x|),itr(|x 3 — x|),...] with “diag” defining a diagonal matrix. 
Let B be the regression matrix with the ith row being hi, where vector 
hi = [1, (Xi - x)^, vn{(xi - x) (g) (xi - x)}"^,...]. 

For any given order P, the parameter /3o provides the estimated function 
value m(X) for a given observation X. We introduce a column vector ei with 
its first element equal to one and the rest to zero so as to extract /3o from 
the least squares minimization problem (Equation]^: 

N 

= ^0 = ef (B^WB)-iB^WY = ^ 1{X)Y,. (8) 

The matrix (B^WB)-iB^W is the Moore-Penrose inverse and ^(X) is the 
so called equivalent kernel of local regression analysis. The parameter P 
determines which polynomial order is used to locally approximate m(X): 
with P = 0, local constant regression is revealed, with P = 1 we obtain local 
linear regression and with P > 1 we obtain polynomials of increasing order. 
Local constant regression assumes the following form for to(X): 

N 

_ Ki\xi-x\)Y, N 

f] it:(|xj-x|) i=i 

i=i 

Equation coincides with the definition of the virial stress in Equation [l] 
Firstly, by choosing tjj to he a uniform kernel function. And secondly, by 
setting Yi to a component of cr^ divided by = V/N (with V being the 
volume of the entire system occupied by Af atoms) forming a stress per atom 
value. In our experimental validation (Section]^ we restrict our attention to 
the LCR (P = 0) and LLR (P = 1) models. 

Remark 1 The equivalent kernel due to LLR may take negative values. Most 
previous related work m in the literature use a non-negative kernel func¬ 
tion. Other references Hills] state the required properties of the kernel func¬ 
tion without explicitly discussing if its either positive or negative. However, 
the kernel function may be allowed to take negative values in general HliEQl. 

Remark 2 The derivation of Hardy stress requires the kernel function to be 
solely dependent on a single distance from the continuum point at x to an 
atom under consideration at x^: ipi = tpiri) with = x^ — x. This allows 
for the desired property d^tpiri) = —dr-tpiri) [55]. Given this, an imple¬ 
mentation of LLR with respect to Hardy stress is not straightforward. The 
equivalent kernel is dependent on both and the distances to the neighbors 
Tj] li{ri,rj). To circumvent this problem, the assumption can be made that 
only the distance is subject of change, if x changes. This can be written as 
9xZ(ri,rj) fj). Therefore, only the distance to atom 

i may change, if the position of x changes. The distance to the remaining 
atoms cannot change. This is one reason, in addition to significantly higher 
computational costs |28j . as to why we focus on the virial stress and do not 
elaborate on the investigation of Hardy stress. 
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4 Numerical examples 

We evaluate the local linear regression (LLR) model on two classical exam¬ 
ples of elasticity. Our first example is the distribution of residual stress about 
the core of an edge dislocation in an elastic plane solid. Secondly, we consider 
an infinite plate with a circular notch under uniaxial load. In both numerical 
experiments we investigate a single crystal of copper in a face centred cubic 
lattice using LAMMPS [11]. This model has 1000 x 1000 x 3 unit cells (with 
length a = 3.615A) with periodic boundaries on the two square surfaces and 
free boundaries on the remaining four surfaces. An EAM potential is used 
with a cut-off radius of 2.5a |31j . Subsequently, energy minimization at zero 
temperature with relative tolerance 10“^® guides the system to equilibrium. 
To model plane strain, the atomistic system is reduced to two dimensions in 
the a;, y-plane by taking a unit cell slice of thickness a in the dimension bear¬ 
ing the periodic boundaries. Furthermore, we constrain our space averaging 
volume to have a circular cross section in the a;,?/-plane. 


4.1 Edge dislocation 


In this experiment we take the base system and create a (100) edge dislocation 
at the centre of the simulation box before energy minimization; this results 
in a Burgers vector oriented along the cc-direction. A circular slice is removed 
from the atomistic system surrounding the origin of the coordinate system 
with a radius of 20o in the x, j/-plane and a thickness of a. This yields 5,024 
data points, the spatial location of which can be seen in Figure 



Fig. 1 (a) The arrangement of the atoms in an edge dislocation gives a rather 
regular sample distribution. Lines of atoms in both x- ((g)) and y-direction (©) are 
selected. A Gaussian (G) and a uniform (U) kernel function for LCR (C) and LLR 
(L) are directly placed above the crossing atom of these two lines. This atom is the 
centre of the shown coordinate system. Horizontal (b) and vertical (c) slices of the 
kernel functions are given for the selected atoms. 


To ascertain the effect of local constant regression (P = 0 in Equation]^ 
and local linear regression (P = 1 in Equation on our atomistic dataset 
we compute the difference between the kernel weights arising from both re¬ 
gression models. For ease of exposition we select a Gaussian (with a standard 
deviation of 4.3lA) and a uniform kernel (with a bandwidth of 7.50A to ad¬ 
just for the different shape of the kernel) [14] for smoothing. The kernels are 
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placed over a selected atoni as shown in Figure and the associated weights 
are computed (slices of the kernel weights are plotted in Figure]^. As can be 
observed, the kernel weights arising from both regression models are effec¬ 
tively identical. The reason for this is due to the almost uniform spacing of 
the atoms which has the effect of rendering the off diagonal blocks in B^WB 
of Equation with values very close to zero in the LLR model. Therefore, 
the kernels for P = 0 and P = 1 are practically identical (see [551 Section 
II.D] for a discussion). 

We take the definition of virial stress from Equation and replace the 
uniform kernel with the LLR equivalent kernel according to Equation [^ and 
set (Ti —>• cTifVi. With respect to the analytical solution of anisotropic linear 
elasticity jS] (elastic constants taken from m), the LLR model obtains a 
0.126 percent increase in root mean squared error (RMSE) as compared to 
the LCR model of virial stress. This difference is negligible and indicates 
that there is no substantial difference between the predicted stress arising 
from models. Replacing the uniform kernel with a Gaussian kernel, thereby 
yielding an altered virial stress, results in a similar observation: the difference 
between the output of both models is barely perceptible (relative increase in 
RMSE: 0.056 percent). 


4.2 Infinite plate with circular notch 


In our second simulation we consider the problem of an infinite plate with a 
circular notch subject to uniaxial loading. The base system is altered by firstly 
removing a hole of radius 60a from the centre following by the application 
of a uniform tensile traction oi p = IGPa on the edges normal to the y- 
direction. The equilibrium state is found at a temperature of zero Kelvin by 
energy minimization with relative tolerance 10“^®. An atomistic solution to 
this problem was previously given in [ni55] . while an analytical solution of 
elasticity theory for the case of plane strain was provided in |13j . However, 
the effect of surface free energy is neglected in the latter as it is considered 
to be small compared to the bulk energy in continuum mechanics. In Section 
4.1 based on a direct comparison of kernel weights, we found that local linear 
regression (LLR) is effectively identical to local constant regression (LCR) 
when applied to an atomistic system that conforms to a regular grid. In this 
experiment we compare both models in the situation where we have a pore 
in the material. This pore introduces a boundary in the system only one side 
of which contains atoms. 

We construct the simulation as follows: a circular slice with thickness a 
centred at a; = 60a and y = OA is taken as our atomistic dataset. The slice 
has radius 30a giving 6,329 data points (Figurej^. We fit both the LCR and 
LLR models to this dataset. The assigned weights for both LCR and LLR are 
shown in Figure]^ for a uniform and a Gaussian kernel (as described in Sec¬ 
tion [T^. If we take the Gaussian kernel as our example, we can immediately 
observe the effect of the data-adaptive equivalent kernel at the boundary of 
the material (Figurej^b)). As there are no atoms in the negative x-direction 
the kernels are therefore truncated in their extent in this region. The red line 
with circular markers denotes the Gaussian kernel weights in the LCR model, 




while the black line with asterisk-style markers denotes the equivalent kernel 
in the LLR model. It is clear that the equivalent kernel has adapted to the 
local density of the data at the boundary: versus the LCR Gaussian kernel, 
the equivalent kernel distributes substantially more weight to the cluster of 
atoms lying on or very close to the boundary (< 2.5A), while considerably 
downweighting those atoms located further away (> 2.5 A). The equivalent 
kernel therefore compensates for the lack of atoms in the negative x-direction 
by up-weighting the contribution of the boundary atoms to the estimation 
of the average stress at the boundary. As the stress at the boundary atoms 
is substantially different to that in the interior, this adaptation is beneficial 
for prediction of stress at the boundary - this estimate will be dominated by 
the stress of boundary atoms, with very little contribution from more distant 
atoms. The net result of this kernel adaptation is a better fit to the atomistic 
data in the boundary region. 



Fig. 2 (a) The arrangement of the atoms at the notch. Lines of atoms in both 
X- (-I-) and y-direction (x) are selected. A Gaussian (G) and a uniform (U) kernel 
function for LCR (G) and LLR (L) are directly placed above the crossing atom of 
these two lines. This atom is the centre of the shown coordinate system. Horizontal 
(b) and vertical (c) slices of the kernel functions are given for the selected atoms. 


Unfortunately, despite our expectations to the contrary, the improved fit 
at the boundary does not translate into an improved estimate of virial stress 
as given by Equation The central cause of this are surface effects which 
have an impact on the atomistic stress per atom values cr^ given in Equa¬ 
tion More specifically, the surface effects cause the stress per atom value 
of the boundary atoms to be very different from the stress per atom values 
of those atoms only a few A into the material from the boundary. This can 
be seen by observing the patterns of the stress per atom values (indicated 
by the diamond symbol) from OA to lOA in Figure |^b). Due to the large 
discontinuity in stress over a small distance in this region, the improved fit 
of the LLR model at the boundary translates into a much poorer fit in the 
region several units of A into the material. The LLR model effectively at¬ 
tempts to smooth out this discontinuity, resulting in a poor fit in the locality 
of the discontinuity. Ideally the impact of the surface-affected atoms should 
be restricted to the surface itself, necessitating that the kernel assign negligi¬ 
ble weight to such atoms when computing the stress at points in the interior. 
However the equivalent kernel, due to its data-adaptive capability, assigns a 
much higher weight to boundary atoms than the LCR kernel, when comput- 
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ing the average stress for atoms in the region adjacent to the boundary (OA 
to loA in Figurej^b)). In contrast, the weights assigned to the surface atoms 
by the LCR kernel do not differ substantially to the weights assigned by that 
kernel to atoms in the interior, effectively curtailing the impact of surface ef¬ 
fects on the predicted average stress of interior atoms. The net result is that 
the classical virial stress (that is, the LCR model), due to its poor fit at the 
boundary, actually yields a superior overall fit to the atomistic data due to 
its insensitivity to the stress discontinuity near the surface. This observation 
negates any potential benefit brought about through modelling the dataset, 
and in particular, the boundary region, using LLR. 



Fig. 3 Stress distribution for notch problem, blormalized normal stresses (a) ctxx 
and (b) ayy in positive ^-direction starting at the radius of the notch and at y =0A. 
The analytical solution (surface effects disregarded) is given by the dotdashed line, 
the volume-related stress per atom values (cri/td) are given by (0). 


5 Discussion and conclusion 

The objective of this paper was to investigate the applicability of automatic 
kernel carpentry for the purposes of modelling atomistic stress using non- 
parametric kernel regression. In many materials of interest, flaws or bound¬ 
aries result in an asymmetry in the distribution of atoms. Automatic kernel 
carpentry provides a principled means of adapting the weighting of the ker¬ 
nel according to the local density of the data samples. It is well known that 
the local linear regression (LLR) model yields a reduction in the bias of the 
regression function, with a modest increase in variance, yielding, in many 
cases, an overall improved fit to the data. This data-adapted kernel is known 
as the equivalent kernel in the statistical literature. In this work we replaced 
the uniform kernel in the virial stress by the equivalent kernel formed by a 
LLR model yielding an altered virial stress. The predictive accuracy of our 
LLR model was then compared to the local constant regression (LCR) model 
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of [28] . The key finding of our experimental validation is that automatic ker¬ 
nel carpentry, as given by fitting an LLR model to atomistic stress, does not 
yield an increase in the accuracy of the virial stress estimate versus the LCR 
model of [55]. The main barriers to effective application of the LLR model 
proved to be both the highly regular arrangement of atoms in the crystalline 
solids under study and the surface-affected atoms at the boundary of the 
materials. 

Firstly, the regular spacing of the atoms in crystalline solids resulted in the 
LLR equivalent kernel to be effectively identical to that obtained using the 
LCR model in the interior. If the atoms are placed in a perfect lattice there 
is no asymmetry in the distribution of the atoms, and therefore no scope for 
kernel adaptation as all points in the space have equivalent density. The net 
effect of this is that the LCR model provides similar stress estimates to the 
LLR model, at a fraction of the computational cost. Secondly surface effects, 
which cause the stress of surface atoms to vary greatly from adjacent atoms 
located in the inner regions of the material, negated any benefit arising from 
the boundary bias correcting capability of the LLR model. There is effectively 
a discontinuity in the stress distribution in the region near and at the surface, 
which contradicts the basic smoothness assumption of non-parametric kernel 
regression. In our experimental validation we observed that the LCR model 
is less effected by this discontinuity than the LLR model. The influence of 
the highly surface-affected atoms at the boundary onto the averaged stress 
in the proximity of the boundary is undesirably magnified by the weight 
distribution of the equivalent kernel, that is, substantially more weight is 
assigned to surface atoms versus those atoms in the interior of the material. 
In contrast, the averaged stress given by the LCR model in the proximity of 
the surface is much less influenced by surface effects. We identified the cause 
of this to be the more equal weighting applied by the LCR kernel to atoms 
in the interior and to those at the surface, thereby providing an average 
that is less sensitive to discontinuities in the atomistic data. In other words, 
the LCR model, due to its increased bias at boundary regions, essentially 
restricts the contribution of surface atoms on the average stress prediction at 
points in the interior. A local linear regression model with change points [23j 
would be a possible means of handling these discontinuities in the atomistic 
data - we leave this investigation to future work. Given these observations 
we found that the additional computational effort of implementing the LLR 
model versus the LCR model is not justified for the modelling of atomistic 
stress in crystalline solids. 

Nevertheless, despite demonstrating that LLR is ineffective for the ex¬ 
ample of crystalline solids, our study has highlighted two useful criteria that 
can be applied by the research community for evaluating when best to apply 
LLR versus LCR. Firstly, if the material under study has a highly irregular 
distribution of atoms in the interior, then LLR is very likely to offer a much 
better fit to the data, and therefore provide more accurate predictions for 
the physical quantity of interest. Secondly, LLR is also likely to give more 
accurate predictions in the situation where the material under study is not 
characterised by wide variations (discontinuities) in the predictive variable 
(in the example of this paper, stress) over small changes in the covariate (in 
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this paper, distance). Given these criteria, we believe LLR may be beneficial 
for similar problems in fluid dynamics and solids with highly irregular struc¬ 
ture (for example, amorphous solids). In addition, the LLR framework may 
also be useful for modelling interfacial regions between two bulk phases (con¬ 
sidered as a two-dimensional continuum) and the transition region near the 
line of contact of three media (considered as a one-dimensional continuum). 
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